scripts/Simlation study/start-simulations.R

In joakimwallmark/PolyOptimalIRT:

source("R/simulateParallel.R")
source("R/scaleThetaFromInf.R")
source("R/scaleThetaToInf.R")
source("R/simulationIteration.R")


# generate using optimal from natmath
load("data/NatMath_fittedmodel.RData")
load("data/NatMath_fittedmodel_mirt.RData") # load natMathGpcm

R <- 1000

# re-estimating models
dataGenerationModels <- list(AnalyzeResult$parList[[10]]$WfdList, AnalyzeResult$parList[[10]]$WfdList,
                             natMathGpcm, natMathGpcm)
dgptype <- c("optimal", "optimal", "gpcm", "gpcm")
thetaValues <- list(AnalyzeResult$parList[[10]]$theta,
                    scaleThetaFromInf(0, 100, qnorm(seq(0, 1, 0.001))),
                    fscores(natMathGpcm, method = "ML")[, "F1"],
                    qnorm(seq(0, 1, 0.001)))
# dataGenerationModels <- list(natMathGpcm, natMathGpcm)
# dgptype <- c("gpcm", "gpcm")
# thetaValues <- list(fscores(natMathGpcm, method = "ML")[, "F1"],
#                     qnorm(seq(0, 1, 0.001)))

cyc <- 15
set.seed(1234)
for (ind in 1:length(dataGenerationModels)) {
    res <- simulateParallel(R, dataGenerationModels[[ind]], thetaValues[[ind]],
                            dataGenerationModelType=dgptype[ind],
                            modelsToCompare=list(c("gpcm", "EAP"), c("gpcm", "ML"), c("optimal", cycles=cyc)),
                            optScr=NatMath_dataList$optList$optScr,
                            usePrecalibrated=F)
    save(res, file = paste("data/simulation-results/n", length(thetaValues[[ind]]), " R", R,
                            " dgp-", dgptype[ind], " precal-F.RData", sep = ""))
}

# pre-calibrated
dataGenerationModels <- list(AnalyzeResult$parList[[10]]$WfdList, natMathGpcm)
dgptype <- c("optimal", "gpcm")
thetaValues <- list(seq(0, 100, 2), scaleThetaToInf(0, 100, seq(0, 100, 2)))
set.seed(1234)
for (ind in 1:length(dataGenerationModels)) {
    res <- simulateParallel(R, dataGenerationModels[[ind]], thetaValues[[ind]],
                            dataGenerationModelType=dgptype[ind],
                            modelsToCompare=list(c("gpcm", "EAP"), c("gpcm", "ML")),
                            optScr=NatMath_dataList$optList$optScr,
                            usePrecalibrated=T,
                            preCalibratedOptimal=AnalyzeResult$parList[[10]]$WfdList,
                            preCalibratedParMod=natMathGpcm)
    save(res, file = paste("data/simulation-results/n", length(thetaValues[[ind]]), " R", R,
                           " dgp-", dgptype[ind], " precal-T.RData", sep = ""))
}


# load("data/simulation-results/errordata.RData")



# simulateResult <- dataSimulation(NatMath_dataList, AnalyzeResult$parList[[10]], nsample=1000)
# save(res, file = paste("data/simulation-results/n", 51, " R", 1000,
#                        " dgp-optimal precal-T.RData", sep = ""))
# simulateResult$mu
# mean(simulateResult$theta[1, ])
# simulateResult$thetastd
# simulateResult$thetabias
# simulateResult$thetaRMSE
# simulateResult$al
# ggplot(data.frame(t=rowMeans(simulateResult$theta), true=seq(0, 100, 2)), aes(x=true, y=t-true)) +
#     theme_bw() +
#     geom_point()
#     scale_y_continuous(limits = c(NA, text), expand = c(0, 0))



# optList <- list(itemLab=NULL, optLab=NULL, optScr=NatMath_dataList$optList$optScr)
# # sample needs to start with 1 instead of 0
# dataList <- make.dataList(NatMath_dataList$U, NULL, optList)
# wfdList <- AnalyzeResult$parList[[10]]$WfdList
# temp <- thetafun(initThetas,
#                  WfdList = wfdList,
#                  dataList = dataList)$theta_out